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Abstract. In this paper, we propose a general way of computing expectation values 
in the zero-range process, using an exact form of the partition function. As an example, 
we provide the fundamental diagram (the flux-density plot) of the asymmetric exclusion 
process corresponding to the zero-range process. We express the partition function for 
the steady state by the Lauricella hypergeometric function, and thereby have two 
exact fundamental diagrams each for the parallel and random sequential update rules. 
Meanwhile, from the viewpoint of equilibrium statistical mechanics, we work within 
the canonical ensemble but the result obtained is certainly in agreement with previous 
works done in the grand canonical ensemble. 



1. Introduction 

The zero-range process (ZRP) is a stochastic model in which many indistinguishable 
particles occupy sites on a lattice [H [2], [3]. Each site of the lattice may contain an 
integer number of particles and the particles hop to the next site with a probability which 
depends on the number of particles at the departure site. In other words, each particle 
interacts with the particles at the same site, i.e., they have a zero-range interaction. 
The hop probability is totally asymmetric, i.e., particles can only move in the definite 
direction and hence the flow of particles, even from the macroscopic viewpoint, never 
vanishes in general. The flow gives rise to a contrasting density distribution of particles, 
and one may therefore regard the ZRP as a typical nonequilibrium system [H E]. 
Since one can choose any function as the hop probability function, the ZRP has been 
extensively studied and applied for a wide variety of many-particle systems PEE]. 

Although there now exist a lot of generalized/extended versions of the ZRP [31 El IS], 
in this paper we focus exclusively on the original one, i.e., the ZRP in one dimension 
and with periodic boundary. Updates of the particle configuration occur at each discrete 
time step, and two typical update rules, the parallel and random sequential update rules, 
are considered. These rules are defined as follows. Parallel update rule: at each time 
step, simultaneously at every site, one particle attempts to hop to the next site with its 
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hop probability. Random sequential update rule: at each time step, a site is selected at 
random and then one particle at the site attempts to hop to the next site with its hop 
probability. 

It is notable that the ZRP in one dimension can be mapped onto an exclusion 
process, i.e., a many-particle system on a periodic lattice whose sites contain a single 
particle or none [3]. Hereafter we call particles in the exclusion process vehicle to prevent 
confusion. According to this mapping, particles in the ZRP are regarded as the distance 
between adjacent vehicles in the exclusion process, i.e., the number of particles at site i 
indicates the number of blank sites in front of vehicle % following vehicle i — 1; accordingly, 
for a particle to hop to the next site in the ZRP is for a vehicle to hop to the next site 
in the exclusion process. This exclusion process adequately simulates a traffic flow in 
which vehicles hop forward each with its own probability depending on the distance to 
the front vehicle. Note that the number of sites in the exclusion process is equal to the 
sum of sites and particles in the ZRP and meanwhile the number of vehicles is equal to 
that of sites in the ZRP. 



2. Nonequilibrium steady state and partition function 

In this section, following [3j, we precisely define the nonequilibrium steady state of 
the ZRP and then formulate the partition function in that steady state. As well as 
equilibrium statistical mechanics, expectation values of physical quantities should be 
given by using that partition function. 



2.1. Nonequilibrium steady state 

Let P({n m }) be the probability of finding the system in a configuration {n m } = 
{rii, ri2, . . • , um}, where n m denotes the number of particles at the mth site. The 
transition probability from {n' m } to {n m }, denoted by T({n m }\{n' m }), is expressed by 
the hop probability function u{n) respectively according to the update rule, i.e., parallel 
and random sequential update rules. Note that it is always true that w(0) = 0. 

We firstly consider the parallel update rule. The transition probability for the 
update rule is given by 
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where, at each site, v m (= 0, 1) particle hops to the next site with probability u(n. 
and S(n) is the Kronecker delta that returns unity if n = and zero otherwise. The 
balance of probability currents at each configuration {n m } is represented by the following 
equation: 



E 

{n' m } 



[T({n m }\{n' m })P({n' m }) - T({n' m }\{n m })P({n m })] = 0. (2) 
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Then, we call a solution of ([2]) the non equilibrium steady state probability (or the steady 
state probability, simply). The steady state probability P({n m }) is given as a product 
of the single-site weights fin): 

1 M 

^(KJ)=~ TT/( ra -) (n 1 + n 2 + "- + n M = N), (3) 

' m=l 

where a normalization Zm,n, the sum of the products over all configurations, is referred 
to as the partition function for the ZRP in the nonequilibrium steady state. One can 
directly confirm that ([3]) is a solution of (121) ; meanwhile the single-site weight is found 
to be 

{l-u(l) (n = 0) 

l-n(l) ft l-ujj) (4) 

Next, we consider the random sequential update rule. In this case, one can describe 
the steady state condition in a simple form: 

M 

2j [u{n m -i + 1)-P( • • • , n m _i + 1, n m - 1, . . .) - «(n m )P({n m })] 0(n m ), (5) 

m=l 

where 6{n m ) is the Heaviside function, which emphasizes that site m must be occupied 
for there to be associated hops out of and into the configuration {n m }. As well as in 
the case of the parallel update rule, the steady state probability P({n m }) is given in the 
product form Q, where the single-site weight for the random sequential update rule is 

{1 (n = 0) 

(See [3] and [10] for details.) 
2.2. Partition function 

It should be noted that f(n) is not identical to the probability that a given site (e.g. 
site 1) contains n particles. Let the probability denoted by p(n), and it is obtained from 
the probability distribution of configurations P({n m }) as 

P( n )= J2 P({n, n 2 , n M }) = fin) M ~ 1 ' N ~ n , (7) 

I 1 AT ^M,N 

"2+"3H hriM=N— n 

The sum of pin) over n is unity by definition, and we thereby obtain the recursion 
formula for the partition functions: 

N 

Z m ,n = Yl f{n)Z M -x,N-n [M > 1, N > 0), (8) 

n=0 

Zi, k =f(k) (k>l). (9) 
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Considering the generating functions, /(C) := g~ =0 f(n)( n and Z M (() ■= Y,n=o Z M,n( n , 
we have a recursion with respect only to M: Z M (() = f(()Z M _i((). Consequently, we 
find the fundamental relation between the two generating functions: 

z M (C)= (hC)) M . (io) 

Thus, the partition function Z M ^ N is obtained from the single-site weight f(n) via their 
generating functions. From the statistical mechanical viewpoint, Zm{() is appropriate 
for us to call the grand- canonical partition function, where we should think of ( as the 
fugacity [H [3] . 

Using p(n) given in <£jty, the average velocity of vehicles in the steady state is defined 

by 

N 

vm,n = S ^u(n)p(n). (11) 

n=0 

In the case of parallel update rule, we find that the average velocity is expressed by the 
partition function as 

EN— 1/ ~\\ri7 

"M,N = -, , N • (1*) 

{—*-) ^M,N 

In the case of random sequential update rule, we also find 

Zm,n-\ /1Q s 
vm,n = • (13) 

(See [12J for details.) 



3. Exact solution 



In this section, we exactly compute the partition function and then the expectation value 
of the average velocity. These computations are fully done by using the formulae for 
the Lauricella hypergeometric function. The definition and formulae used are relegated 
to Appendix. 



3.1. Parallel update rule 

We start with the hop probability function u{n) given as 
u(0) = 0, 

< u{n) < 1 (1 < \fn < K), (14) 

u{n) = l (\/n>K+l), 

where K (> 1) is a constant integer. This is not a strong restriction, for one may take 
K to be as large as necessary. As far as traffic-flow models are concerned, it is quite 
natural to choose this kind of hop probability function. 
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Substitution of (THj) into (J1J yields f{n) — for n > K + 1, and we accordingly 
transform into a convenient form, i.e., 

k+i k r « 1 /.v" 



c-n 



j=1 u(j) 



(15) 



7(c) = £/wc n =(i-u(i))(i+oE 

n=0 n=0 

where we define nj=i cij = f if n = 0. For further computations, we formally factor the 
most right hand side of (TT5"j) : 

K 

/(o=/(o)(i+on( i +^o- as) 

i=i 

Then, comparing (fT5|) with (fl6|) . we note the relation between xi, a?2, ■ ■ ■ , xk and 
u(l), u(2), . . . , u(K), namely 

e (x) = l, e n (x) = TT 1 ~^[ J) (n = l,2,...,K), (17) 

where e n (x) is the nth elementary symmetric function of K variables, x = {x\, . . . , xk)- 
(Note that the elementary symmetric functions are defined by its generating function: 
ELI 1 + x iC) = J2n e n( x )C n , and specially e (x) = 1.) 
Accordingly, substitution of <HM into ( fTUl) yields 

_ftT 

^m(o = /(o) M (i + o m n( i + x & M - 

Compare ffl8|) with flA.21) . the generating function for the Lauricella hypergeometric 
function Fq, and one finds an expression for the partition function in the nonequilibrium 
steady state: if N < KM + M then 

K 

Zm,n = /(0) M (-l) Ar ^^F D (-iV, ^M~^M, M-N + l;x), (19) 



and otherwise Zm,n — 0. (See Appendix A for definition of the notation.) As seen 
below, this expression using the Lauricella hypergeometric function is quite useful for 
exact computation of expectation values. 

In the case that iV > KM + M the average velocity in the steady state takes the 
value of unity, since all vehicles take more than A'-site distance from the front vehicle 
(or all sites in the corresponding ZRP contain more than K particles). If N < KM + M, 
by substituting of ( |T9|) into ( |T2l) we have an expression for the average velocity. Using a 
recursion formula for the Lauricella hypergeometric function with respect to parameters 
dAppendix A.2[ ), 

aF D (a + l, (3, 7+I; x)-^yF D (a, (3, 7; x) = (a-j)F D (a, (3, 7 + I; x), (20) 
the sum in the numerator in (fl2l) is firstly carried out: 
(-M) 

zZ Vy^ F d(-^, -M,...,-M,M-n + l;x) 

= M)jV ~ 1 F D (l - N, -M, . . . , -M, M - N + 1; x). (21) 
UJjv-i 
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Accordingly, the average velocity of vehicles for arbitrary numbers of sites and vehicles 
is given as 

= N F D (1 -X.-U ■ l:.rj 

v m,n M Fd(—N, —M, . . . , —M, M — N + 1; x) 1 ' 

We remark that (1221) is valid for arbitrary numbers of M and N. 

Then, we expand (1221) in lattice size L = M + N with a fixed density of 
vehicles p = M/L. Our strategy to take the thermodynamic limit is to eliminate the 
hypergeometric function from fl22|) . The Lauricella hypergeometric differential equation 
that Y = Fo(a, /3, 7; x) satisfies is defined by 

[( rf -l + 5)5i-x i {a + S){j3 i + 5i)]Y = (i = 1, 2, . . . , K), (23) 
where 5 = J^^t an d $i = Xid/dx{. Using a differentiation formula (1A.11I) . we have 

F D (1 - N, -M, . . . , -M, M - N + 1; x) = y - ^6y, (24) 



where y = Fq(—N, —M, . . . , —M, M — N + 1; x), the denominator in ([22]) . and then 
substitution of this into f[2"2l yields 

N N 
by = -r{h-v MjN )y (h = — ) (25) 

Meanwhile, we consider the hypergeometric differential equation for y (where a = —N, 
Pi = -M and 7 = M - N + 1) and find 

x N VM,NXj + %(1 - Xi)5iV M , N 

Siy = — ■ y (i = l,2,..., K). (26) 

h 1 - V M ,N + V M ,NXi 

Since 5iy = by, equating the sum of the right hand side of (126"]) with that of (125]) we 
successfully eliminate y: 



h r,, v ^ ; ■ (27) 

i=l 

Expand the average velocity as fjw.iv — ^0 + viL^ 1 + v 2 L~ 2 + • • •, and substitution of 
this into (126]) yields a series of equations for v , V\, v 2 , ■ ■ ■■ Consequently, as N tends to 
infinity, we have 

h-v = Y-^ , (28) 

1=1 

where v (= Vq) is the average velocity in the thermodynamic limit. Solving the series 
of equations for v\, v 2 , ■ ■ ■ in turn, one can obtain the higher-order correction terms as 
functions of v. The computations are straightforward, but however the expressions of 
those terms are too long to be included. 

In jXQJ , they already obtained the same result as ( 128]) working in the grand-canonical 
ensemble where the fluctuation of the particle number N is allowed. Actually, (]28]) is 
consistent with their result, 

(N)=w-^logZ M (w), (29) 



/V V MjN Xi + jf(l - Xi)5iV M ,N 
^ 1 - V M,N + V M ,NXi 
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where (N) means the expectation value of the particle number in the grand-canonical 
ensemble. Since the right hand side of (1281) is in a symmetric form with respect to 
x\, . . . , xk, one can express it using the elementary symmetric functions of them. Recall 
Ylifi + x i() — Ylm e n( x )( n , an d one finds the identity, 

wxj = J2r=i re r (x)w r 

The left hand side of (1301) is equivalent to the right hand side of (1281) if one lets 
w — v/(l — v), and thus we find that (1281) is equivalent to (1291) . 

We note that using ( 1291) one has the fundamental diagram in a parametric 
representation: flux Q{y) = p{v)v and density p{v) = 1/(1 + h), where v (0 < v < 1) 
is the parameter. As an example, we illustrate a fundamental diagram choosing a hop 
probability function that approximately describes motion of vehicles in traffic flow, i.e., 

tanh(n — c) + tanhc , N 

u(n) = : , 31) 

v ; 1 + tanhc v ; 

where we let c = 3/2. Note that the hop probability function can be estimated from real 

traffic data, and (13~T1) is often chosen for traffic-flow models [11], [13]. Figure H] shows the 

exact fundamental diagram illustrated by using (12"§1) with K = 50, and the simulation 

data in addition. They show complete agreement. 



0.25 




Density 



Figure 1. A fundamental diagram of the exclusion process corresponding to the zero- 
range process with the parallel update rule. Gray line shows the exact solution given 
in the parametric representation, and black points are simulation data. Numerical 
simulation is done with site L — 1000. These are in good agreement. 



3.2. Random sequential update rule 

In the same way as in the case of the parallel update rule, we obtain the fundamental 
diagram for the ZRP with the random sequential update rule. 
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The single-site weight is given in ([2]), and then the generating function for these 
weights is given as 



n=0 



n 1 



(32) 



where u(n) is given in fTH]) . Then, we can factor (132]) into f(() = (1 — Q 1 Y\f=i{^-~ x iC) 
and express the partition function by the Lauricella hypergeometric function: 



A' 



'M.N 



(M) N 



F D (-iV, -M, . . . , -M, -M - N + 1; x), (33) 
where the independent variables x = (xi, x 2 , ■ ■ ■ , xk) satisfy 

n 1 

e (x) = l, e n (x)=TT — (n = l, 2, ...,#). (34) 

From ( JT3l) and ( 1331 . the average velocity for the finite numbers, M and N, is expressed 
by the Lauricella hypergeometric function: 

= N /pi-.V • 1- -M -M. -M - A • 2: .,•) 

wm.tv y . v - 1 F D (—N, —M, . . . , — M, — M — N +1; x) ' 1 
Thu>. applying the formulae for the hypergeometric functions we find 

h _ N _ f 1_\ V MjN XiVM,N + gjg)g ^V f36 , 

M V MJ1-V M , N ^ l-XiV M ,N 

As well as the parallel update case, we can expand vm,n as a series t>o+f \L~ l +v 2 L~ 2 +- • •. 
Especially in the thermodynamic limit, we recover the result given in[5J E] 

(N)=v^-logZ(v), (37) 

which was obtained within the grand-canonical ensemble. 

In contrast with the parallel-update case, if K tends to infinity, we need to determine 
the range of parameter v. For this purpose, we impose a simple condition on the hop 
function, i.e., 

lim u(n) = (0 < Mqo < 1), (38) 

n— >oo 

which surely allows us to determine the parameter range. Then, the range of values of 
v becomes < v < Moo- It depends on the asymptotic behaviour of u(n) whether f(v) 
converges or not at v — Uoo. The convergence of f{v) is discussed in the context of 
condensation occurring in the ZRP [31 HH [15] . 

We choose a traffic-flow model with the hop probability function (13T|) and random 
sequential update rule, and show the exact fundamental diagram and simulation result 
in figure [2j These are in good agreement as well as in the case of the parallel update 
rule. 

We consider a simple hop probability function which allows us to provide a direct 
formula of fundamental diagram instead of the parametric representation, i.e., 

u(l)=A, u{n)=p (n>2). (39) 
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Figure 2. A fundamental diagram of the exclusion process corresponding to the zero- 
range process with the random sequential update rule. (See Fig. Q] for comparison.) 
Gray line shows the exact solution given in the parametric representation, and black 
points are simulation data. Numerical simulation is done with site L — 1000. These 
are in good agreement. 




Substitution of PJJ into fl37]) via ([32]) yields 
h 



p 2 v 



(40) 



(p ~ v)(Xp + (p — X)v 

Accordingly, eliminating v we obtain the fundamental diagram explicitly described by 



Q = pp 



A(l-q)p(l-p) 
2(l-g)(l-p) 



where q 



X 

P 



(41) 



One can also figure out an explicit formula for this model in the case of the parallel 
update rule; however, one has to solve a cubic equation. 

If we take the limit as A — > p (i.e., q — > 1) in (141)) . we recover the fundamental 
diagram of the ASEP, Q = pp(l — p), which is a well-known solution [16]. If we let 
p = 1 (q = X, accordingly), then we recover another exact solution given in [17] , where 
they obtained the same result by using the matrix-product ansatz. 



4. Summary and remark 

In this paper, we provide an exact solution of the zero-range process: we have 
the partition function for the general case, which is expressed by the Lauricella 
hypergeometric function, and which allows us to make further computations. Although 
the solution obtained is rather formal, that is enough for us to compute expectation 
values of macroscopic quantities. In particular, using the exact partition function we 
obtain the average velocity of vehicles for arbitrary numbers of sites and particles. 
Moreover, expanding the average velocity in terms of the site number, we obtain the 
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fundamental diagram as the leading term of the expansion. One can also figure out the 
correction terms through a long computation. 

Due to the fact that the steady state of the ZRP has a factorized form, one 
can do within equilibrium statistical mechanics while dealing with the macroscopic 
quantities averaged over all the lattice. In this regard, one may take Zm,n as the 
canonical partition function in which the particle number is conserved. As a result, we 
compute the expectation value of the average velocity within the canonical ensemble. 
In general, statistical treatment in the canonical ensemble is more difficult than that 
in the grand canonical ensemble. Both results in these ensembles are equivalent in the 
thermodynamic limit. An analytical approach based on the grand canonical ensemble 
was previously made to the ZRP with the random sequential update rule in [6], and 
to the ZRP with the parallel update rule in |10j . Recently, canonical analysis of the 
condensation is also done [HJ[T5]. In this regard, it will be interesting to exactly compute 
the probability p{n) that a single site is to contain n particles in the steady state. 

Recently, as a traffic-flow model, we introduce the stochastic optimal velocity (SOV) 
model [TT|[T8] in the form of combining two exactly solvable models, i.e., the asymmetric 
simple exclusion process (ASEP) [16J appearing as a special case of the ZRP, and the 
ZRP itself; however it does not seem to us that the model has an exact solution in the 
general case. In numerical simulations, we find that the SOV model has an intriguing 
and complex phase transition in the fundamental diagram, and our future goal is to 
analyse in detail the phase transition observed, e.g., by using a theoretical method for 
exactly solvable models. 
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Appendix A. The Lauricella hypergeometric function 

There are four types of the Lauricella hypergeometric function, being denoted by Fa, 
F&, Fq and Fq. Each of them has a series representation, an integral representation and 
a differential equation it satisfies. Since only Fq appears in this paper, we confine our 
attention to F D . (See [19] for details.) 

Firstly, we define the Lauricella hypergeometric function Fq with K arguments by 



where a, Pi, ... , [3k and 7 are all complex parameters, (a) n = a(a + 1) ■ ■ ■ (a + n — 1) 
is the Pochhammer symbol, and notations in the right hand side are defined as follows: 
Ik = {tti = (m%, . . . , m K ) rrii E Z> (i = 1, 2, ... , K)}, then, for m G Ik, 




(A.l) 
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\m\ = m 1 + \-m K , (P) m = ((3i) mi ((3 2 ) m2 ■ ■ ■ (/3 K ) mK , (l) m = (l) mi (l) m2 • • • (l) mK 

and x m = Xi mi ■ ■ ■ xx mK ■ We also use the following notations for simplicity: (3 = 
(Pi, Pi, Pk) and x = (x x , x 2 , . . . , x K ). 

The hypergeometric differential equation for F D (a, P, 7; x) is given in (j23|) . We 
omit the integral representation because it is not used in the present work. 



Appendix A. 1. Generating function 

We present the generating function for Fq: 

K 00 / \ 

(1 - C) K II^ 1 - = E [ -^F D (-n, P, K — n + 1; x)C- (A.2) 

i=l n=0 

Here, we give a proof of flA.2j) by mathematical induction. One should note the relations, 

00 ( — \ 

n=0 

and 



q(o+l)---(o + Z-l) (-l) m (a)/ r . 

( a + /-m) ••> + /-!) (l-a-Z)m' 1 1 



(i) If if = 1, then one has 



(i-cr(i-*iO-* = E E 7, wff "^ 



m ™ /-k+m\ 



mi=0 fc=0 

00 00 



(l)fc(l)mi 



= E -T^M-n, Pi, k - n + 1, sOC* 1 - (A.5) 

n=0 ^ ' n 

Thus, ([Oil holds if if = 1. 
(ii) Assume that (1A.2j) holds if if = I — 1, and one accordingly has 

(i - cr n( 1 - *<cr A = (i - *icr A E K^M-n, p'^-n+v, *'k 



(i) 

1=1 n=0 v ' 



-ro^-'/t, p , n — 11 -|- j.; x 

n 



yy^ r^L fc+n (-n)\ m \(p ! )m(Pl)k(x') m Xl h 

,n n W» 77" («-n+ l) H (l) m (l)* 

fc=0 n=U mHi_i 

^ W« («-» + i)N + *(i)-(i)* ' 

where x' = . . . , P' = (Pi, . . . , Pi-i) and I k = {m = (mi, . . . , m fc ) ; m 8 G 

Z> (i = l, 2, fc)}. 

Thus, ( 1A.2I) holds also when K = I. Therefore, we may conclude that (1A.2j) holds for 
any natural number K. 
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Appendix A. 2. Recursion formula 

We provide a recursion with respect to the parameters, a and 7, for F&: 

aF D (a + l, P, 7+1; x) = 7-F D (a, j3, 7; x) + (a-~/)F D (a, P, 7 + I; 2). (A. 7) 

One should note simple formulae for the Pochhammer symbol, 

n 1 1 ft 
(a+ l) n - (a) n = -(a) n , and — — — — = — . (A. 8) 

a (a) n (a + l) n (a)„,+i 



From (lA.ip and the two formulae above, one respectively has 

F D (a + 1, p, 7; x) - F D (a, (3, 7; *) = £ ^ ^M^ ^ ( A .9) 

_ OL CY) | m | ^1 J m 



me/ 

and 



F D (a, /3, 7; *) - F D (a, /5, 7 + 1; x) = £ ~ rXTT^ ( A ' 10 ) 



7 (7+l)|m|(l), 



Then, we can construct (I A. 71) from flA.91) and (jA.lOp . 



Appendix A . 3. Differentiation formula 

Since the Lauricella hypergeometric function is a multivariable function, when taking 
the derivative one needs to deal with partial differentiation: 5i = Xid/dxi and 5 = J2i 
In the present work, we specially use the following formulae: 

5F D (a, 0, 7; x) = a (F D (a + 1, P, 7; x) - F D (a, P, 7; x)) 

= (7 - 1) (F D (a, P, 7 - 1; x) - F D (a, p, 7; x)) . (A.ll) 

Note that 5iXi mi = rriiXi" 11 , and from (1A.8j) one finds 



5F D {a, p, 7; x) = 



m 1 



\m\x 



X 



a ((a + 1) H - («)h) (/?) 
fel (7)M(l)m 

rf , 7 ~! - r^ 1 ) {a) ii p)m x™ (A.12) 

^Wy-Dh ( 7 )h/ (l)m 



Thus, (lA.lip is proven. 
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